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Abstract 

The problem of determining the neuronal current inside the brain from measure- 
ments of the induced magnetic field outside the head is discussed under the assump- 
tion that the space occupied by the brain is approximately spherical. By inverting 
the Geselowitz equation, the part of the current which can be reconstructed from 
the measurements is precisely determined. This actually consists of only certain 
moments of one of the two functions specifying the tangential part of the current. 
The other function specifying the tangential part of the current as well as the radial 
part of the current are completely arbitrary. However, it is also shown that with 
the assumption of energy minimization, the current can be reconstructed uniquely. 
A numerical implementation of this unique reconstruction is also presented. 

1 Introduction 

Magnetoencephalography (MEG) is a non invasive technique that can be used to 
investigate brain activity. The physiological basis of MEG is the following: The main 
functional units of the brain are certain highly specialized cells called neurons. For higher 
mental processes the most important part of the brain is its outermost layer called cerebral 
cortex, which contains at least 10^° neurons. When neurons are active they produce small 
currents whose basis is the change in the concentration of certain ions p[j (ionic currents). 
The flow of current in the neural system produces a weak magnetic field. The measurement 
of this field outside the brain and the estimation of the current density distribution that 
produced this field is called MEG. Other names often used are magnetic source imaging, 
magnetic field tomography, and current-flow imaging. 
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Neuromagnetic signals are typically 50-500/T, which are of the order of 10~^ of the 
earth's geomagnetic field. Currently, the only detector that can measure these tiny fields 
is the superconducting quantum interference device (SQUID). The theory and practice 
of SQUID as applied to MEG measurements, as well as several practical approaches for 
shielding all other external magnetic fields except that of the brain, are discussed in the 
excellent review f2]. Here we only note that the SQUID, which is the most sensitive 
detector of any kind available to scientists jHj, is based on the exploitation of several 
quantum-mechanical effects, including superconductivity as well as the Josephson effect. 
The SQUID can be thought of as a digital magnetometer where each "digit" represents one 
flux quantum, and it is essentially a transducer converting a tiny change in magnetic flux 
into a voltage. Whole-head magnetometer systems are now used by several laboratories 
in Europe, USA and Japan. 

The current density J and the magnetic field B are related by the Maxwell equations. 
These equations can be simplified using two facts. First, the permeability of the tissue 
in the head denoted by /i is that of the free space, i.e. /i = fiQ. Second, the quasistatic 
approximation is valid, namely the terms dE/dt and dB/dt can be neglected, where E 
denotes the electric field and B denotes the magnetic induction^. Using these facts the 
Maxwell equations become 

V-B = 0, VAB = /ioJ, (1.1) 

where • and A denote the scalar and vector product, respectively, and V denotes the usual 
gradient. Part of J is due to neuronal activity, and part of J is due to the electric field E, 

J = JP + ctE, (1.2) 

where Jp denotes the neuronal current (primary current) and a denotes the conductivity. 
The electric field E satisfies V A E = 0, thus there exists a scalar function V, called the 
voltage potential, such that 

E = -VV. (1.3) 

Making the further assumption that a = aj inside the head and a = ao outside the head, 
where ao and aj are constants, equations ()1.H) - ()1.3|1 imply the celebrated Geselowitz 
equation |lj 

- ^{ai-ao)f y(y)n(y) A x ^ fi, 

where |x| denotes the length of the vector x, Vl denotes the volume occupied by the head, 
dVL is the boundary of fi, n denotes the unit outward vector normal to the surface dVl, 

^Let (7 and e denote conductivity and permitivity which are assumed to be uniform, and let E = 
Eo(x) exp(27ri/t), where / denotes frequency. Then Maxwell equations imply that the term dE/dt can 
be neglected provided that \edEi/dt\ <C ItrEj, or 2Tr fe/a <C 1. This is indeed the case since for the brain 
a = 0.3r2~^m~^, e — lO'^eoj ^^nd since in neuromagnetism one usually deals with frequencies of about 
lOOHz m, 2TTfe/a ^2x 10"^. Similar arguments hold true for the B field. 
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and dS denotes the infinitesimal surface element on d^l. For a recent rigorous derivation 
of this equation see jH]. 

Equation ()1.4j) relates Jp inside the head with B outside the head. However, it also 
involves the value of V on the surface of the head. This serious complication can be 
avoided if one makes the simplifying assumption that the head is spherical. Then, and if 
in addition o-q = 0, which is justified since o"o -C c/, equation ()1.4|) reduces to [HI- [El 

B = /ioVf/, 

1 / . J-WAy^y (1.5) 

47ry|y|<i |x-y|(|x||x-y| +x- (x-y)) 

Equation ()1.5|) relates Jp inside the head (|x| < 1) with B outside the head. This equation 
is the starting point of many of the algorithms used in MEG. It defines the following inverse 
problem: Given B, which is obtained from the measurements, find Jp. 

The main difficulty with the above inverse problem is that it is not unique. This fact 
was already known to Helmholtz since 1853 9J. For example, it is clear from equation ()1.5p 
that the radial part of does not contribute to U. However, in spite of intense scrutiny 
by many investigators, the fundamental question of which part of Jp can be reconstructed 
remained open. 

Here we first give a complete answer to this question, see theorem 2: Jp can be 
uniquely decomposed in the form 



1 fdG 1 dF\ 1 r 1 dG dF . 
p \o9 smO dip J p \sm9 dip o9 ' 



where e^, eg, are the unit vectors associated with the spherical coordinates {p,6,(f) 
and J'^, G, F are scalar functions of (p, 6, (f). This decomposition for vector fields on the 
sphere is the analogue of the celebrated Helmholtz decomposition for vector fields on R^. 
We will show that knowledge of U determines only certain moments of F with respect to 
p, while and G are arbitrary. More precisely, it can be shown that U can be represented 
in the form U = J2eni'^^,rnP~^^^^^y£,miO,(f), where Yi^m are the usual spherical harmonics 
and the constants are determined from the measurements. Then we will show that 
F can be represented in the form F = ^ fi,m{,p)Yt,m{6, V'), where only the moments of 
fi^m are determined in terms of q^^. 



p'^'fi,Mdp=i'2i+l)ce,m. 

The above results imply that by decomposing Jp into a "silent" component and into 
an "effective" component, we can show that the Geselowitz integral operator provides an 
one to one map of the effective component of Jp into the magnetic field B, or into the 
magnetic potential U, outside the brain. Furthermore, given U the effective component 
can be explicitly computed. We emphasise that, since the decomposition into a silent and 
into an effective part is of a general nature independent of any assumptions on Jp, our 
result that U determines the effective component of the current uniquely and says nothing 
about the silent component, is actually a general statement which is model independent. 
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The next part of the paper deals with the case when we assume some relations between 
the effective and the silent components: We will show that if one requires that Jp is such 
that energy is minimized, then Jp is indeed unique, see theorem 3: In this case J^, G, F 
are given by the equations 

r = G = ^, ^ = EE (1.7) 

1=1 m=-e 

In addition to the above analytical results we also present a numerical algorithm which 
given U (p, 9, ip) for one specific value of p > 1 and for some equally spaced values {Oi}}^'^'' 
and {</?j}o""''') it first computes and then computes Jp using equations ()1.6|) and ()1.7|) . 

The non uniqueness of the inverse problem has been the Achilles heel of MEG. For 
example in the most comprehensive review on MEG 2 , it is written that "with the 
assumption that MEG mainly reflects the activity in the tangential part of the cortical 
currents" , while in ^0! it is written "what cannot be seen should not be looked for" . Even 
the "father" of MEG, D. Cohen has stated JT] "identifying those tangential sources, rather 
than localization, is the real use of the MEG, there is no localization magic". We hope 
that both the analytical and the numerical results presented here will contribute towards 
determining the advantages as well as the limitations of MEG. 

Regarding other brain imaging techniques we note that at present time the most 
important such techniques are the functional magnetic resonance imaging (fMRI) the 
positron emission tomography (PET) and the single photon emission computed tomog- 
raphy (SPECT), as well as a new version of electroencephalography (EEG). These tech- 
niques involve tradeoffs among the following important considerations: temporal resolu- 
tion, spatial resolution, invasiveness, and cost. Assuming that the question of uniqueness 
of the MEG is answered, the spatial resolution of MEG (1cm), of PET and SPECT (4- 
5mm), and of fMRI (1.5mm) are similar; the spatial resolution of the conventional EEG 
is quite poor. On the other hand the time resolution of EEG and MEG is much better 
than that of PET, SPECT and fMRI. The time resolution of PET, SPECT and MRI is of 
the order of 1 second, while that of MEG and EEG is of the order of 10 milliseconds. This 
is a crucial factor if one wants to study brain dynamics. For example, MEG data suggest 
that speech areas of the brain are activated 100 milliseconds after the visual areas. MEG 
is the only truly non invasive method. EEG is minimally invasive (placing electrodes on 
the scalp), while in PET, SPECT and MRI the subject is exposed to radioactive tracers 
and to strong magnetic fields, respectively. EEG requires a rather inexpensive apparatus 
(of the order of thousands of dollars). The fMRI has the advantage that can be obtained 
by modifying the existing MRI apparatus. PET employs positron-emitting radionuclides 
which have such short half-lives that there must be a cyclotron near the site of scanning, 
thus the cost is of the order of multimillion dollars. The cost of the MEG is similar to 
that of the PET. 

We conclude the introduction with some remarks: 

(a) We expect that the combination of our analysis of the spherical model and pertur- 
bation theory can be used to study realistic head geometries. In this respect we also 
note that progress has been recently made regarding ellipsoidal geometry [T^ . 
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(b) The question of what additional information can one obtain by measuring E (using 
EEG) is under investigation. 

(c) Due to the orthogonahty of the decomposition of into silent and effective com- 
ponents, the assumptions that the norm of the solution is minimal, implies that 
the silent component vanishes. Clearly, one can assume other relations between the 
silent and the effective components, for example one may assume that the current 
consists of a finite number of dipoles. It is well known that this assumption, un- 
der certain conditions, also leads to a unique solution. This current can also be 
represented in the form (jl.fjj) with F of a particular form, and therefore can be con- 
sidered within our formulation. Thus the answer becomes model dependent only at 
the stage when one makes an assumption about the form of Jp. For other models 
see [lH!-[in]. 

(d) "Least-square" methods have been used extensively in inverse problems. However, 
our approach of using such methods in order to find an approximate numerical 
solution of the Geselowitz equation is fundamentally different than the existing 
ones. Indeed, it is based on the explicit decomposition of the current into a silent 
and an effective component, and thus could not have been used before obtaining 
this decomposition. 

(e) In practice, the magnetic field is measured approximately over a half-sphere over the 
head and not over a whole sphere. However, since in the numerical reconstruction we 
assume a finite number of spherical harmonics, the approximate knowledge of U over 
part of a sphere is sufficient to determine approximately the current. Clearly the 
problem becomes more and more ill-posed when the number of spherical harmonics 
increases. A stability result in this direction is under investigation. 

(f) It has been correctly pointed out by one of the referees that it is sufficient for 
the solution of the inverse problem to invert dU / d\x\ instead of U. Furthermore it 
has been correctly pointed out that this latter inversion is much simpler since the 
expression for dU /d\^\ is simpler than the expression of U (see 

(g) A short summary of the analytical results presented here was announced in 

2 Analytical Results 

We first show that equations ()1.5|) can be written in an alternative form, which is more 
convenient for determining the part of which can be reconstructed from the knowledge 
of f/(x). 

Theorem 1. Let f/(x) be defined in terms of Jp by equation ()1.5|) . Then t/(x) can also 
be expressed by the alternative representation 

^(x) - ^ £{(V. A JP(z)) . zU .^IzHzl) .y, |x| > 1. 

(2.1) 
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Proof. Let /(z) denote the following function of z 



m = n / ' ^ ^ ^) ■ (Vx$(x))}_M.c^|x|, (2.2) 

|Z| Jo 1^1 

where $(x) G C^(]R^). We will integrate J(z) over the sphere |z| < 1: We first multiply 
by |zp and integrate with respect to c?|z| along < |z| < 1. Interchanging in the resulting 
expression the order of the integration with respect to d\x\ and to (i|z| we find 



1 ri / rl 

2 



/(z)|z|^ci|z| =47r {Vx$(x)}^^M, ■ ^ ^ JP(z) A |x|z (i|z| J (i|x|. 

We then integrate this equation with respect to rfz, z = z/|z|, and denote |x|z by y. This 
yields 

/ /(z)rfz = -47r f $(y) Vy ■ f ^ f ^y) A y|z| d\z\\ dy. (2.3) 

^iz|<i j|y|<i \|y| J\y\ viy| / / 

It is straightforward to show that 
Indeed, the rhs of this equation equals 

where cp denotes cyclic permutation; the Ihs of equation ()2.4j) equals 

d 

oy3 

and using the chain rule as well as noting that several of the resulting terms cancel we 
find the expression ()2.5|) . 

Using equation ()2.4|) . as well as noting that the term Vy(|y|~^) is perpendicular to 
JP A y, the rhs of equation ()2.3|1 becomes 



-4n I <|.(y)^ ( A(V. A JP(z)) ■ z}^_^ Iz^lzA dy. 
J\y\<i \y\ \J\y\ / 

Replacing the rhs of equation ()2.3j) by this expression and replacing /(z) by the definition 
()2.2|1 . equation ()2.3j) and the standard Green's function representation for solutions of 
Poisson's equation, give equation 1)2.11) provided that the result of the lemma proven in 
the appendix A is valid. Note that according to our proof, equation ()2.ip is valid in the 
distributional sense, but simple regularity arguments imply that it is also valid pointwise. 
QED 
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Theorem 2 (Representation theorem) 
The vector J^(x) can be uniquely decomposed in the form 

JP(x) = JP{p, e, ^)e, + J\p, e, ip)ee + J'^(p, 6, ^)e 



(2.6) 



where ep, e^, e^ are the unit vectors associated with the spherical coordinates p > 0, 
0<9<7i,0<Lp< 2tt, and the scalar functions and J'^ can be represented in the 
form 



1 fdG 



1 dF 



1 



1 dG OF 

+ 



p \ sin 6 dip 89 



p \d6 sin 9 dip 

where G(p, 6', ip) and -F(p, 6', yj) are scalar functions of the arguments included. 
Assume that f/(x) is defined in terms of by equation ()1.5p . Then 



(2.7) 



1 



f 1 


f ^ f 


^|y|<i 








y 


2 /i 



^ / Ae,^F{\z\,e,^)d\z\]dy, 
M -'lyl / 



|x| > 1, 



(2.8) 



where Ag^^ denotes the Laplacian with respect to the spherical coordinates 9 and ip, i.e. 



A, 



sin 6' 



d 



39 



d 



— sin 6'— + 



89 



sin 9 dif"^ 



Proof. We first decompose Jp into a radial and a tangential component. Clearly gives 
no contribution to U. Also the tangential component can be uniquely decomposed in the 
form ()2.7|) . see appendix B. Using equations ()2.6|) and ()2.7p we find 

1/19. dF 1 d''F\ 

|z| \sm9d9 o9 sm 9 Oip^ J 
and ()2.H) becomes equation ()2.8p . QED 
Corollary (Non uniqueness of the inverse problem) 

Assume that U{x.) is defined in terms of by equation (jl.Sj) . Let a vector Jp(x) be 
written in the form ()2.(ij) where the scalar functions and J'^ are given in terms of the 
scalar function G and F by equation ()2.7j) . 

The function f/(x) is independent of and of G, and furthermore only certain mo- 
ments of F can be computed in terms of U. In particular, F{p, 9, ip) is given by the 
expression 

oo i 

F{p, 9,^)=Y,Y1 krn{p)Y,,UO, V), P < 1, < < TT, < ^ < 271 , 

i=l m=-i 

where Y^^m are the usual spherical harmonics, the moments of fe^m{p) can be determined 
in terms of q^^, 

i [ p'^'kMdp = (2^ + l)ce,m, (2.9) 

^0 
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and the constants can be determined from the given data using the fact that U (x) 
can be expressed in the form 

oo e 

U{p, 9,^) = Y,Y. ce,n.p-'^'^'%m{0, P > 1, 0<9<n, < ^ < 27r. (2.10) 

£=1 m=-£ 



Proof. Equation ()2.8p imphes 



AU 



Ae,^F(|z|,^,^)rf|z|, |x| < 1, 



|x| 



(2.11) 



AU = 0, |x| > 1. 
Let us represent F and U in terms of spherical harmonics by 

F{p,9,Lp) = '^fi,mip)Yi^miO,ip) and U{p,9,Lp) = ^ n<?,m(p)>£,m(6', v?)- 



i,m 



Then equations ()2.11|) imply 



2 , ^(£ + 1) 

p p^ 







p > 1, 



where prime denotes differentiation with respect to p. The general solution of the homo- 
geneous problem is ap^ + Pp'^^'^^^ where a and P are constants and £ is a positive integer. 
Since m as p oo it follows that 



-(^+1) 



To solve the inhomogeneous problem we use variation of parameters in the form U£^m{p) 
Ai^rn{p)p^- This implies 



(4,™P''+')' = -^(^ + 1)«£,™(P), «£,™(P) = / / kmip'W- 

Thus 

A[^^p''^' = i{i + 1) I ae,UpW + A'e,^a). 
Convergence at p = implies 

A',^^{1) + i{i + I) I ae,Mdp' = 0. 
Using Ai^ra = Ui^mP~\ wc find 



(2.12) 



p=l 



i2£+l)c£,m. 
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This equation together with ()2.12|) imply 



+ 1) / aLMdp= (2£+1)q,, 



Using integration by parts we find ()2.9|) . 

Theorem 3 (Minimization of energy) 

Define the energy by 

W = 



QED 



(2.13) 



x|<l 



Then if 

JP = rep + J'ee + J'^e^, 

where and J'^ are given by equations ()2.7|1 . it follows that the minimum of W under 
the constrain 



£=1 m=-£ 



'£,ra{p)dp = {2i+l)ce. 



where Y^^^ are the usual spherical harmonics and are given constants, is achieved 
when 



(2.14) 



=1 m=-£ 



Proof. Substituting equation ()2.(jp and ()2.7|1 in the rhs of equation ()2.13|1 we find 



W 



x|<l 



1 fdG 



dG 



+ 



1 f dPy 1 fdF\ 



- 



\d6 J p2 gj^s Q yg^p J p2 gj^2 Q yg^p J \d6 J 



- 



(ix. 



where we have used that the term involving G^Fq — GgF^ vanishes, 

•1 r-T r2n 1 



JO 



sin 6 



dG dF dG OF 

86 dtp dtp do 



p sin 6dpd6dip = 0. 



The constraint involves only F, thus it follows that the minimal energy is achieved 
when JP = G = and when H is minimal, where 



H 



1 /"TT /'27r 

Jo Jo 



1 fdF 



p2 sin^ 9 \dip J p^ \dO 



1 fdF 



p^ sin OdpdOdif. 



(2.15) 



The term inside the bracket equals |VFp — (§^)^, which using integration by parts (with 



either F or ^ equal to at |x| = 1), equals — [FAF + (ff )^], where 



^ ^ d^F 2dF 1 ^ 
op'^ p op 
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Using 



£,mi 



and the orthogonality of the spherical harmonics, it follows that 



H -- 
Hence, 



£,Tn 



f'lmip) + -firnip) 
P 



fe,m{p) 



kmip) + ifupyr p'dp- 



Thus, provided that either /^,m(l) or /^'^(l) equals zero^, we find 

H = Y,K£ + l) ['fUp)dp. 

The assumption that /^,m(l) = is without loss of generality since the tangential part of 
the energy which is given by equation ()2.15p does not involve differentiation over p, thus 
in general ()2.15|) can be obtained by approximating / by functions equal to zero at p = 1 
and then passing to the limit. 

The minimization of this H, under the constraint ()2.9j) . implies ()2.14|1 . 

We note that equation implies that f/(x) behaves like 0(p"^), hence £ > in 

equation ()2.10p . cqo = 0, and the sum ()2.14|) starts with £ = 1. QED 

3 Numerical Implementation 

In equation ()2.1()j) „ denotes the spherical harmonics, namely 

= a,,^P,,^(cos^)e™'^, 



> 1, < m < 



where the bar denotes complex conjugate and 



I2£ + 1 {£~m)\ 



Pe^m are the Legendre functions, namely 



(3.1) 



(3.2) 



Jm 



where 



Pf(x) 



1 



ix' - lY 



2V. dx^ 

are the usual Legendre polynomials of degree £. 

For the numerical implementation we replace in the sums appearing in ()2.10|) . 
oo by £max, where £max is chosen by the procedure explained below. 

^These conditions are true since the support of JP lies in the interior of the sphere. 



(EH 
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3.1 Computation of q 



TO 



We first discuss how to compute ce^m from either U{p, 9, ip) or from B(p, 9, ip). 
Suppose we know U{p, 9, ip) for one specific value of p > 1 and for some equally spaced 
values 9i, (pj, such us 

<9i <7r, i = 0,.. .,imax, 
<iPj < 2TT, j = 0,. . .,jmax- 

Using the orthogonality of Y^^m equation ()2.1()|1 implies 

Ui9,ip)Ye,miO,^)dicos9) dip = ct,mP'^'^^^ ■ 



-1 / 

Therefore, using the first equation in fj3.1|) . we obtain 

r27r / /'TT 



ij U {9, ip)Pi,micos 9) sin 9d9] e-'^^'^dip. 
Using (jSUl, we find 



Ci,-m = i-irc^, i>l, 0<m<i, 



(3.3) 



where 



and 



/»27r i'2'K 

Ue^m = / Ui^mi^) cosmifdif — i / tj^^m{ip) sinnnpdip, (3.4) 
Jo Jo ' 

Ui,mi^)= / f/(^^,^)P£,„^(cos^)sin^rf^. (3.5) 
Jo 

For the numerical calculation of the three integrals appearing in ()3.4|) and ()3.5|) we 
use an extended closed formula, namely 

f{x)dx = Ax Q/i + + ^/s + /4 + • • • + fn-3 + + + ^/n j • 



For the numerical calculation of the Legendre functions P^^m(cos^) we use subroutine 
plgndr from Numerical Recipes [T7]. The constants ai^m are given by ()3.2|) . 

Suppose we know B = {Bi, B2, B3) instead of U. Then, using the first relation in p.5j) 
and spherical coordinates we obtain 

Up = sin 9 cos ipBi + sin 9 sin (pB2 + cos 6' S3, (3.6) 

where 

~^^HpA^^ . = 1,2,3. 
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Moreover, by differentiating ()2.10|) with respect to p we find 

U,{p, 9,^) = -J2Yl + ^)P~^'^'%UO, (3.7) 



=1 m=-e 



Thus, if we know B, we can compute Up from (jH.fij) and then we can compute Ci^m from 
()3.7|1 . following the same procedure as before. 

The choice of Imax- Using ()3.H) and the second relation in ()3.3|) . the real part of 
(|TTn|l implies 



Ir, 



f/(p,^,^) = (3.8) 
i=i 

Re{ce,o)ae,oPe,oicos6) + 2 ^(i?e(Q,m) cosmip - Im{ci^rn) smm(p)ai^rnPe,micos9) | . 



m=l 

Differentiation with respect to p yields 



f/,(p,e,^) = -£(£+l)p-(^+^)- (3.9) 
-Re(Q,o)a£,o-P£,o(cos6') + 2 ^(i?e(Q,m) cosm^? - Im{ce^rn) smmip)ae^mPe,m{cos9) 



m=l 



Therefore, after calculating the coefficients q^^, following the procedure outlined earlier 
we can use either ()3.8p or ()3.9|) to re-evaluate either U or Up. In this way not only can 
we test the efficiency of our procedure, but we can also run our program several times, in 
order to find the most appropriate value for imax- 

3.2 Computation of the Minimizing Current 

Using relations ()2.14|) . ()3.1|) and the second relation in ()3.3|) . the real parts of the 
functions J^, J"^ defined in fl2.7|) (with G = 0) are given by 

A.,«,.)^-^E^?^ii^/. (3.10) 

m{Re{cE,m) sin rmp + Im{ci^rn) cosmip)ai^rnPe,rnicosi 

\m=l 

and 

n,,e,,,.-.,ne'f:^^^±^M±^/. (3.11) 

Re{cefi)ae,oPloicos9) + 2 ^(i?e(Q,m) cos my? - Im{ci^rn) sinm(/})a^,„P^' ,„(cos6') | . 

m=l / 
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Recall that the Legendre functions satisfy the recurrence relation 



Therefore 

-sin ePI,Q{cos9) 
-sin^P(„(cos^) 

Thus, in order to calculate numerically the current we apply the following procedure: 
We take some 6 and if points, such that 0<^<7r,0<(y9<27r. We first calculate 
the Legendre functions Pi^m{cos9). In a separate subroutine we calculate the quantities 
Pe,m{cos9)/ sin 9 (for this purpose we have developed a subroutine similar to plgndr). 
These quantities appear in both and the second relation of ()3.12j) . Note that these 

quantities are valid even for ^ = or = tt. We then calculate from ()3.12|) the quantities 
— sin 6'P/„^ (cos 9). Finally, we take a value of p such as < p < 1 and calculate J^(p, 9, ip) 
from ()3.10p and J'^{p,9,ip) from ()3.1H) . In all the above we use the imax value that was 
found with the procedure outlined in the previous subsection. 



P£,i(cos 
mcos9 
sin 9 



for m = 0, 
■Pi,micos9) + Pi^rn+i{cos9), for m > 0. 



(3.12) 



3.3 Verification of the Algorithm 

We have tested our numerical algorithm for several functions U{p,9,(p). In what 
follows we discuss two typical examples. 

Example 1 

Let U be given by 

U (p, 9, (f) = —2 cos 9— + sin 9 cos 9 cos ip— — sin^ 9 cos 2o9— . 

p2 pi pi 

Note that this function has the form ()2.1()|1 with = for £ > 2. 

First, we evaluate U ioi p = 1.5 and some equally spaced 9i and ipj, where imax = 100, 
imax = 200. We calculate numerically the coefficients from the first relation of ()3.3|) . 
()3.4|1 and ()3.5j) . and then evaluate from ()3.8j) the approximate value of f/, at the above 
p, 9i and (pj. Furthermore, we start with Up instead of a f/, we calculate in a similar 
way and then calculate the approximate value of Up from ()3.9|1 . 

We run our program several times with Imax from 1 up to 40 and we found that the 
best value is Imax = 2, which is consistent with the exact form of U . For this value the 
difference \U — Ua\ is of order 10~^, at most. 

Secondly, we calculate numerically J^, J"^ , using ()3.1()|1 - ()3.12|1 . in the above 9i, ipj and 
some equally spaced p^, such as < p^ < 1, namely A; = 0, . . . , kmax, where kmax = 25. 
Then we calculate analytically c^^m from ()2.10|) . F from ()2.14|) . and J^, J*^ from ()2.7|) . For 
the above U we have 

35 35 
F = — 30p^ sin 9 + —p^ sin 6* cos 6* cos — —p^ sin^ 9 cos 2^9. (3.13) 
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The analytical and the numerical values of J and J"^ in the various ^j, ^pj and pk are 
almost the same (the absolute value of their difference is of order 10~^, at most). 

We have also verified the validity of equation ()2.8j) as follows: We take F from (jH.lHj) 
and evaluate numerically Ua from ()2.8|1 : \U — Ua\ is of order 10"^, at most. 

Example 2 

Let U be given by 

^ 1 pixi + P2X2 + P3(a^3 - a) i^N 
47r + + (X3 - a)2]3/2 ^"^^ ^ 

with a = 0.5 and (^1,^2,^3) = (0.1, —0.2,0.6). We evaluate U for p = 1.5 and the same 
equally spaced 6i and cpj, as in Example 1. We again calculate numerically the coefficients 
Ci^rn and then Ua- 

For this example we found that the best value for Imax is 10. For this value the 
difference \U — Ua\ is of order 10~^, at most. 

Finally, in Figure 1, we present the density plots of the minimizing current ( J^)^ + ( J*^)^ 
for the above function U in various cuts perpendicular to the xa-axis. 




Figure 1: Density plots for the minimizing current of the function U given by ()3.14|) . 
Starting from top left = -0.9, -0.8, -0.6, -0.4, -0.2, 0, 0.2, 0.4, 0.5, 0.6, 0.8 and 0.9. 
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Appendix A 

Lemma. Let 

Ui^, = f ^"j . X, (A.l) 

^ ^ V|x-z|(|x||x-z| +x- (x-z))y ^ ^ 

Then 

/ (A,t/(x,z))$(x)dx = /'''{(J(z) Ax) ■ (Vx*(x))K=|xj z (i|x|, (A.2) 

jRf |Z| Jo 

where A is the Laplacian (i.e. A = V ■ V), and $(x) G C^f (M^). 

Remark. As AU is singular close to x = the integral in the Ihs of ()A.2jl should be 

understood in the sense of distributions. 

Proof. Let z be at distance a from the origin along the direction Xg. Let l^e(z) denote a 
small neighborhood of the interval [0, z] defined as follows, 

a(z) = C,(z) U5,(0) U5,(z), 

where C^{z) is the cylindrical region 

Ce{z) = |x' e : p= \Jx[^ + x'2 = e, < 4 < a| , 



while 5*^(0) and ^'^(z) are the semi spherical regions 

^,(0) = {x' G : |x'| = e, Xg < O} , 

and 

5,(z) = {x' G : |x' - z| = e, Xg > a} , 

respectively. 

Let $(x) be a test function, then from the theory of distributions it follows that 
At/($) = / (At/(x,z))$(x)cix = / f/(x,z)A$(x)rfx = 

lim /" ?7(x, z) A$(x)cix = - lim /" f f/^ - rf^, (A.3) 

where dS" denotes the infinitesimal surface element on the surface dQ^{z), n denotes the 



unit outward normal, and we have used the fact that AU = in R /(^^(z). Let /i(z,e), 
/2(z, e), hi^z, e) denote the contributions from the integration along Ce(z), S^{0), and S'e(z), 
respectively. It is easy to show that lirn/2 = lirn J3 = 0. We now compute Ii. Let 

/(x', z) = |x' — z|(|x'||x' — z| + x' ■ (x' — z)). 

Thus if x' G Ce(z), 



f=[{a- X',f + + ^/ 2 ^ + _ ^/ )2(^2 ^ ^/ 2 _ 
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Hence 



p 



p2 + X- 



+ 2pVp2 + (a _ a;/^)2 + 



+ _ ^/ )2 



and 



dp 



2 2Vp^ + 4^ 



P2 + 4 



v/p2 + (a - x'J^ 



+ Pf, 



where / is bounded at p = 0. Evaluating /, and at p = we find 



ap2 



/ 



p=0 



4(« - 4)^ + (a - 4) (4^ - «4) = 0' 



9p2 



5/ 

dp 

/ , («-4)^ 



p=0 



2a;; + 



p=0 



+ 2{a - x'^) + 



/ 2 _ I 2 



3 a X3 0.3 

i9p dp 

x' = (4>4)4)> z = (0,0,a), J = (Ji, J2, J3), 



(J(z) A z) • x' = a{J2x[ — x^Ji) = ap{J2Cosip' — Jisimp') 



The integral ()A.3|) involves —U^ + ^W-- Also, since 



it follows that 



where we have used x'l = pcos(f' and = p sirup'. 

Equations ()A.1|) and ()A.3|) imply that we need to compute 



lim 



2tt p\z\=a 




pd(p'dx'^[a{J2 COS (f' — Ji sin (/)')] < —-7^ — h $ ( -7 — 



^0 



p (9$ 



/5p 



1 p a/ 



However, $(x') = $(pcosv9',psin(y9',X3), thus as p — > 0, 



(9$ 5$ 
$ = <l>(0,0,x;j) +p cosy.'— (0,0, x;,) +p sin y.'— (0,0, x;,) +0(p2). 



9x'j^ 



and 



(9$ (9$ 9$ 



(A.4) 



(A.5) 



Substituting the expressions for $ and for |^ in ()A.5|) . it follows that the rhs of equation 
[T3I) involves 



'•271- ^^3 



ap'^ df 
p^^^Jo P dp 



lim 



,9$ 

9x']^ 



9$ 

5X2 



9$ 



p3 a/ 



dx'i J P^O /2 c?p ' 
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But 



Also 



Thus, 



where we have used 



pdf 



hm 

P-o / dp 



hm 



dp ^ ^ap2 



9p 



1 + hm 



ap2 P dp^ 

ay 
ep2 



hm ^ = hm ^ 



p' df 



2x', 



hm -r- = 4 hm — = 4Ht 



P-O /2 5p 
Hence 



pp 



hm/i 



47r 
a 



9$ 



9$ 

^2(z)^(0,0,x:,) 



In the above derivation we have used the convenient set of coordinates x', such that 
z is along x'^. This result can be immediately generalized by writing Ji in an invariant 
form. Then (IA.2I) follows. 



Appendix B 

We will show that and J'^ can be expressed by equation ()2.7|) . Indeed, if 

J = J'ee + J^e^, 
then the corresponding 1-form on the sphere of radius p is 

p psmO 

On a compact Riemannian manifold, any 1-form a has the unique decomposition 

a = dG+ (-1) *(i + 

where G is a function, /? is a 2-form, a'^ is a harmonic 1-form, and * is the Hodge operator. 
Also there do not exist any nonzero harmonic 1-forms on the sphere. Furthermore, 
*(3 = F, where F is a function. Hence 

a = dG+i-l)*dF. 

Using 

, „ dG dG 
dG = —d9 + —dip, 
oO OLf 

and 

1 OF OF 
*dF = -T—^^de - smO—dif, 
sm 6 dip oO 
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we find 

09 sin 9 d(f' d^p 09 ' 

and equations ()2.7|) follow. 

Remark. In the case of M^, the analogous decomposition is given by Helmholtz 
theorem: Let A = A^i + A^j + A^k, where i, j, k are the unit vectors along the x, 
y, z axis, be a vector field in M.^. Then there exists a function G and a vector field 
B = B^i + B^j + B^k such that A = VG + V A B. A relationship between the general 
decomposition and the one in M.^ can be established using the following facts: (i) A 
differential 1-form a = dx + dy + dz can be canonically identified with the vector 
field A, where A^ = a^, = a^, A^ = a^. (ii) In the Hodge operator transforms 
a differential 1-form a into the differential 2-form (3 = (3^^ dxdy + (3^^ dydz + (3^^ dxdz, 
where (3^^ = a^, f3^^ = —a^, (3^^ = a^. (iii) There do not exist any nonzero harmonic 
1-forms in M^. 
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